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The random-phase approximation with second-order screened exchange (RPA+SOSEX) is a model 
of electron correlation energy with two caveats: its accuracy depends on an arbitrary choice of mean 
field, and it scales as 0(n 5 ) operations and C(n 3 ) memory for « electrons. We derive a new algorithm 
that reduces its scaling to 0(n 3 ) operations and 0(n 2 ) memory through controlled approximations 
and derive a self-consistent mean field that approximates Brueckner coupled-cluster doubles (BCCD) 
theory with RPA+SOSEX, referred to as Brueckner RPA (BRPA). The algorithm similarly reduces 
the scaling of second-order M0ller-Plesset (MP2) perturbation theory with smaller cost prefactors 
than RPA+SOSEX. We test the new model's accuracy on the dissociation of H2. 



The random-phase approximation (RPA) combined with 
second-order exchange models the correlation energy of the 
uniform electron gas in the high-density limit 1 . For densities 
relevant to materials, this model is extended to second-order 
screened exchange (SOSEX) within coupled-cluster theory 2 . 
RPA+SOSEX matches quantum Monte Carlo benchmarks on 
the uniform electron gas 3 to within 0.002 Ha/electron. Further 
extension to finite and inhomogeneous systems is possible by 
combining the RPA+SOSEX model with Kohn-Sham orbitals 
and orbital energies within density functional theory 4 (DFT). 
Benchmarks on finite systems 5 show mean absolute errors of 
0.012 Ha/atom for total energies of first and second row atoms 
and 0.008 Ha/molecule for atomization energies of molecules 
in the G2-1 test set. Unfortunately, B3LYP-based DFT is more 
accurate than RPA+SOSEX, even after further corrections^. It 
is unclear whether the dominant source of error is the use of a 
Kohn-Sham mean field or the RPA+SOSEX model itself. The 
optimal RPA+SOSEX mean field is as yet unknown. 

Despite a growing interest^ in RPA methods, including 
RPA+SOSEX, DFT remains the computational workhorse of 
electronic structure theory because of the high cost and poor 
scaling of RPA calculations. DFT requires C(n 3 ) operations 
and 0(n 2 ) memory for n electrons, while 0(n 5 ) operations 
and C(n 3 ) memory are required by RPA+SOSEX. Even for 
small n, RPA needs a larger basis than DFT for convergence 
of sums over virtual orbitals. When n is large, RPA+SOSEX 
is computationally intractable with existing algorithms. 

In this paper, we propose a new RPA+SOSEX algorithm 
to match the computational scaling of DFT and define a self- 
consistent mean field. The algorithm uses a spatially localized 
basis, fast summation of the Coulomb kernel, and a low-rank 
approximation of energy denominators-^. The RPA+SOSEX 
model is again revised to approximate a higher level of theory: 
Brueckner coupled-cluster doubles (BCCD) theory 1 1 . BCCD 
applies a nonunitary similarity transformation that produces a 
non-Hermitian mean field. We replace it with a best Hermitian 
approximation to preserve the structure of orthogonal orbitals 
with real orbital energies that is expected of a mean field. 



We consider a Hermitian many-electron Hamiltonian with 
n electrons over an orbitals in second quantization notation, 

H = E + h xy a.\ d y + \ V xy n x h y -\V xx h x , ( 1 ) 

for raising (d\) and number (n x = d\d x ) operators in a basis 
with indices {w,x,y,z} and V vv = V yx . Repeated indices in a 
tensor are summed implicitly unless they appear unrepeated 
in a tensor in the same equation. We define a reference Slater 
determinant, |$), with an orbital transformation, 

d x = <j) px c p , (f)* px (j) qx = 5P, c}\$) = c a \$)=0. (2) 

General orbital indices are {p,q,r,s}, occupied orbital indices 
are {i,j,k,l}, and virtual orbital indices are {a,b,c,d}. With 
Sqx = ^px^qxi the H coefficients in the orbital basis are 

^ = </>;aAv, v q p r = spy XY sf y . (3) 

We antisymmetrize tensors with Vq S r = Vq" —Vqi '. 
BCCD theory approximates the ground state as 

|*)«rap(f 2 )|*), f 2 = \T t f ctc4cj. (4) 

We enforce partial symmetry constraints: Tj" b = Tj] a but not 
T t f = -Tff. %f is defined by projection of H\V) =E\^!), 

($|Xexp(-f 2 )//|*) =£($|iexp(-f 2 )|*), (5) 

forX e {l,cl c a ,c jCbCi Co}- These BCCD equations are 

E-E =\Qt i +l$ = li+\yi + \v£f!j', (6) 

0— U a 7? a b — \7ak.Hrcb .nrac^jkb ,rfac\iklrpdb 
~~ ®i ' K ij ~ V ic J kj +I ik V cj +I ik V cd 1 lj ■ 

0^f+b:T i f-b*T$+b b c Tr-b%?+R<* 

. 1 Trabnncd 1 1 nrabxjkl 1 1 Hnabxjklnncd 

+ 2 V cd 1 ij + 2 1 kl v ij + 4 J kl V cd 1 ij i 

where b p q is a non-Hermitian Brueckner mean fielcP-Q Eq. (|6]) 
needs 0(a 4 n 6 ) operations and 0{a 2 n 4 ) memory to solveP. 
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We derive RPA+SOSEX by approximating T t f with 

T f „ ^ V ;-\ s- = 5?° + SgV^frjad, (7) 

for w a j = e fl - £j > 0, and retaining all terms in Eq. |6]i that 
preserve this form. Only the fourth-order tensor equations are 
modified, to an antisymmetrized RPA Riccati equatiorfSl, 



= V i f + (u Jai + uj b j)T i f+D^ 

r\ab _ -y/ak^pcb .rj^ac-y/kb .rj^ac-y/kl rj^db 
U ij ~ V ic 1 kj +1 ik V cj +1 ik V cd 1 lj ! 



(8) 



which is unsymmetrized by removing ' ' to define T t j fully. 
The rest of Eq. Q is an extended Hartree-Fock (HE) theory, 

E = E a +\p xy {h yx +fy X ), f xy (j>py = e p (j) px , (9) 

= K K y + V x 5 X y ~ P X yV X y + S^, , 

Pxy = 4>ix4>*iy , Pxy = 5 xy - p xy , V x = V xy p yy , 

that includes a Hermitian correlation correction term, 

£ „ = \ (CTxy + CT* X + PxzVzwPwy + Pxz^tzPwy) , (10) 

a X y = 4> ax cj>i{y b +v xz sf z -sf z v zy )f l f. 

We refer to this model as Brueckner RPA (BRPA). 
Eq. (jTji enables Eq. ([8]l to be factored about V xv , 

SW^jy = (S^+T^S^WxyiS^ + T^S^y). (11) 

S" x is defined by equating the two factors and expanding 7y , 



Sl = Sf x -Sly r S b Xj(u: ai + Lo bi ) 



cancels from this equation to define A xy (oj ai ), 

A%{u J ) = -S<SS%/{u Jai + u J ), 
A„.(uj) = B xy (uj)+A x: (uj)V m .B W y(uj), 

B xy (u) = A%(u) - SfV^AUadSl/iUai + W), 

which has an analytic continuation from uj ai to uj G C. 
We split (uj ai + Lo bj )~ l inA° (w ai ) and B xy (oj ai ) using 



(12) 



(13) 



X xv (uj ai ) = 



dfi x xy (n) 

2iti uj„i — 



forXe{A°,fi}, (14) 



cm i cin 1 G yx (ri)G xv (n') 



2ni J r 27r/ 51 — f2' + ko 



with closed counterclockwise contours, r v separating e a from 
e, - icj, and r o separating e, from e a + iuj, and the mean-field 
Green's function, G xy (uj) = 4> px 4>* Y /(ui - e p ). We calculate 
these integrals approximately using numerical quadratures, 

l/(u> ai + LL> bj ) W tt e /[(uj ai + Ul e )(uJhj-UJ e )], (15) 

for /?i points, w e , indexed by {e,/,g}, /?2 points, u> a , indexed 
by {a,b}, and p\ points, w/, indexed by {/,.;'}■ A' an d fh 
have logarithmic error dependence and weak n dependence^. 
We combine a and / into one index, p. We use a permutation 
of indices, e, to manifest symmetry: = u* and = f2 e . 



We approximate Eq. ( p"3j ) by applying a closure relation, 

l/Kwoi+w^^ai+a;/)] « /((j^+Ug), (16) 
Af = -(^-5/)/ (Wc - W/ ), Af =-V>, 

where V* are coefficients of a finite difference approximation 
of ^(waz + cj)" 1 for oj = uj e . This induces approximants, 

X xy (cu ai ) » (l e Xt/(u ai -u.) forX G {A°,A,B}, (17) 



when combined with Eq. ( 15 1. The approximant coefficients 
are AJ° = A° (w e )> B xv = B xy(Ue), and A* that solve 



A^-Q^G^, Gjy = G xy (ujp), (18) 
A%, « lA,l', ; /}; iv - fi/A/A/ V W B^, 
pa A* + 7 Af A^KwAf, , 



zx z w wy i 



which is a reduced form of Eq. ( p"3j ) for 



Eqs. (1 15b and 



16 1 reveal reduced numerical rank in T- 



ab . 



(19) 



Tab ^ f~<a fL r\e jrefr\f fb_ „J_ 
1 ij ~ ~ U ax U xi il qi U xy i>j ' l-v 'v; • 

u% = c:ln g v zw c{,z , C# = 5}5 xy + n g Af/v xz Aj y ■ 

Tensor hypercontractiorPSuses a similar form without orbital 
energy dependence: Tf- b w -</>*^, x t/ xv 



We apply Eqs. ( 12 1, ( 15 1, and ( 19 1 to simplify <t xv to 



°xy ~ ^ .„<'., - ' -i^z - ^xw^wy ^wz , 

= f2 e V vz A_ lv ,y vvv , y^; = G- xz h zw Gzi,y, 

k% = ^cx+^nifii^f^iy, 



(20) 



4 = ny3%w^+nip f M G§ z u e Js§ x 
- [Yjj'ij-ni;,/'' + nip^u^i. 

The first (second) term in contributes the RPA (SOSEX) 
correlation energy to E in Eq. (|9|l through S AV . The analogue 
in second-order M0ller-Plesset (MP2) perturbation theory is 

ifiPfL _ r ia jfaP — o r 1 * C- v 



vy xy ^ ' xz zw wy ? J xv — a '« i4 ^ u ig' t W Lr M5i 
= + ^qp^z^zy- 



'xzy ' 



(21) 



m xy yx ai xz zy 1- 'zyx 

-^yy yz j^+n: A G%y zy jfy, 



combined with HF orbitals and energies defined by replacing 
fxy4> P y = £ P 4>px in Eq. ^ with (/ u >:».,)• '.» = e p (j) px . 

Like HF and DFT, BRPA is a self-consistent field theory, 
with the added burden of calculating a xy . In a conventional 



algorithm, we iteratively solve Eq. ( 12 1 for S" x and calculate 



<7 AV from Eq. (10 1. In a fast algorithm, we iteratively solve 



Eq. (18 i forAl, and calculate <j„ from Eq. (20i. We assign 
Px i-> VxyPy a cos t of a~fn operations, where 7 is reduced by 
a fast summation method. Both algorithms improve upon the 
0(a 3 n 5 ) operations and (D(a 2 n 4 ) memory cost of the known 
RPA+SOSEX algorithm^. Table ^summarizes their costfP. 
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TABLE I. Leading-order, per-iteration floating-point costs of tensor 
formation 16 for two MP2 algorithms and two BRPA algorithms. 



Tensor Operations 


Memory 


Conventional MP2 


2aV 


a% 4aV 


2 2 
an 


Conventional BRPA 


2a 2 n 2 


S? 4ctV 


2 3 
on 


<r rv 6ct 3 r? 


2 2 


Fast MP2 


2a 2 rc 2 + a(A>+/3 3 )« 


A% c^Aftftrc 2 


ia 2 /3m 2 


J% 2a 3 /3 2 ftn 3 + 2a 2 /3 I /3 2 ft» 2 


a 2 /3 in 2 


4° 2a 3 /?! (2ft + 2ft + 7 )« 3 + 8a 2 ft l3 2 /3 3 n 2 


a 2 (f3 2 + l3 3 )n 2 


a% 2a\(3 2 + p 3 )n 3 


a 2 n 2 


Fast BRPA 


2aV + a(/3,/3 2 + A/93)n 


A'„. 6a 3 ftn 3 +2a 2 /3i(3A+7)n 2 


a 2 /3,« 2 


B' xy 6a 3 / 9 1 « 3 +2a 2 /3 1 (3/3 1 +7)n 2 


a 2 /3,n 2 


U$ 0(a 3 /3 2 n 3 + a 2 /3 2 (/3,+ 7 )n 2 ) 




7| 2a 3 /3 2 ft« 3 




4, 2q 3 [(4ft + 7 )/3 2 /3 3 + 0\(.I32 + ft)]n 3 


Q 2 (ft + ft)« 2 


cr,, 2a 3 (/3 2 + ft)» 3 


2 2 



The conventional and fast algorithms also apply to MP2. 
Intermediate calculation of cr® enables orbital optimizations 
such as in the BCC2 model 17 . By utilizing fast summation of 
V xy , the fast MP2 algorithm improves upon the 0(a 4 n 4 ) cost 
of the previous best algorithm^. Additional approximations 
such as orbital localization 19 and stochastic s ampling^lSU can 
reduce MP2 costs even further and achieve linear scaling. 

We test BRPA on a small problem that is difficult for both 
DFT 22 and RPApl the dissociation of H2. H2 is an extended 
Hubbard model (a = n = 2) in the zero-differential-overlap 
(ZDO) approximation 24 , which fits the form of Eq. dill for 



/' 





-t 


" 




u 


u 


V 


V 





1" 





-t 


, v = 


u 


u 


V 


V 


-1 





A* 







V 


V 


u 


u 


_ 


— / 





M _ 




_ V 


V 


u 


u _ 



Its HF, MP2, BRPA, and exact (BCCD) total energies artP 

E HP = E + 2(ji-t) + ^U + V), (23) 
p _p (u-vf p _ F (u-vf 

E = Eo + 2fi+4\(U + V)-^4t 2 +j\(U-V) 2 . 

With total energies of HJ, E + = Eq + /i - 1 , and H^, we 
invert {E+,E ,Ehp,Empi} to determine {/i,t,U,V}. These pa- 
rameters are fit to accurate energies^ to evaluate the BRPA 
solution semiempirically in Fig. [T] All approximations fail at 
dissociation. Significant reduction of the orbital energy gaps 
improves BRPA, consistent with smaller RPA+SOSEX errors 
in calculations based on semilocal DFT 7 compared with 
The BRPA mean field approximates quasiparticle gaps^S that 
are systematically underestimated by semilocal DFT. 

BRPA is a poor approximation of BCCD theory. A change 
of mean field can reduce errors in the total energy, but further 
reduces fidelity to the underlying theory. A variational theory 
of RPA+SOSEX optimized over mean fields may resolve this 
discrepancy, but first it must repair known instabilities 27 . 




2 4 6 8 10 12 14 1.0 1.5 2.0 2.5 



R (Bohr) R (Bohr) 

FIG. 1. Total energy, E, and virtual-occupied orbital energy gap, uj, 
of H 2 as a function of interatomic distance, R. Shown are the exact 
(black), HF (red), MP2 (orange), B3LYP (green), and BRPA (blue) 
results. For R < 4, RPA+SOSEX is exact with a tuned gap (cyan). 
Gaps are model-derived, with the exact gap from BCCD theory, and 
a noninteracting gap (magenta), 2f , indicative of semilocal DFT. 

Efficient BRPA calculations must minimize a, (3\, Pi, ft, 
and 7. Minimizing ft, ft, and ft is a well-defined problem 
of optimizing integration quadrature. Minimizing 7 requires 
further development of fast Poisson solvers. To relate Eq. ([T]) 
to the standard quantum chemistry Hamiltonian, we interpret 
x as a real space index and use an a — > 00 limit to recover the 
continuum. <p px is approximated in a finite basis that restricts 
orbitals in number and accuracy. BRPA will converge slowly 
in the size of this basis, similar to other electron correlation 
models. There are multiple interpretations and solutions: as a 
slow virtual orbital summation it is corrected with remainder 
modeling 28 , and as a cusp in the many-electron wavefunction 
it is corrected with cusp modeling^!. With Gj v calculated as 
shifted matrix inverses of f xv and p vv calculated with 

ft y = / ^Gf v (fi)wnK4» (24) 
Jr„ 27r; - - 

the direct use of orbitals, including virtual orbital summation, 
is removed from BRPA. Basis convergence generalizes to an 

operator approximation problem for G% and Ai,. Similar to 

p 

wavefunction cusps, Gj y and AL both have singularities that 
are poorly approximated in an orbital basis. They need direct 
modeling to reduce the basis set size required to approximate 
the remaining smooth operator. Direct approximation of Gj y 
and Ai, also enables use of operator-specific basis sets with 
more flexibility than a single atomic orbital basis set. 

Fast RPA calculations enable a compromise between DFT 
efficiency and coupled-cluster accuracy. The C(n 3 ) cost of an 
RPA correlation model is much greater than the 0{ri) cost of 
a standard DFT model and much less than the 0(n 6 ) cost of 
the BCCD model. It matches the (D(n y ) scaling of mean field 
calculations that generate the physical inputs needed by these 
models. Whether an RPA model can surpass the accuracy of 
standard DFT models and approach BCCD accuracy remains 
to be determined. Several open problems also remain before 
efficient implementation of RPA calculations on real materials 
attain practically useful performance. 

I thank Jay Sau, Norm Tubman, Jeff Hammond, Andrew 
Baczewski, Rick Muller, and Toby Jacobson for many useful 
discussions. Sandia National Laboratories is a multi-program 
laboratory managed and operated by Sandia Corporation, a 
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SUPPLEMENTARY MATERIAL 
I. TENSOR FORMATION 

Efficient evaluation of tensor equations requires a careful 
ordering of operations and choice of intermediate variables. 
To account for the costs in Table [l]of the paper, we reduce the 
calculations to pseudocode where each addition, subtraction, 
multiplication, and division is a floating-point operation and 
each floating-point number is a unit of memory. X and Y are 
workspace tensors that are reindexed as needed. Workspace 
sizes are listed next to the algorithm names in Table [I] 

Conventional MP2: 

<-0 
for each a do 



jx 



' V xy<t>jy<l>% 



X 

for each i do 

Yjx *~ Xjxfiix 



JX 



Yjx 

j • x f/(e a - 

Y, <- Yj x p jx 
X^h[X) 

X x <- V xy Y y 
o , _o i j, 



ei + e b -e:) 



a. 

end for 
end for 

Conventional BRPA: 
for each a and i do 

Xx «- (old) 
-Z^(old) 



b* y (X+X X -Xy) 



xl 
xl 



■ Xf/(e a - ei + e b -e j) 



X 



JX 



b.x 



^(new). o:,o n \X n o] x 
end for. 



(T xy <- 

for each a do 

Xjx^VxyS a jy 

for each i do 

Xj ^~ Xj X S^ x — Xj X S b j X 
Xf^Xf/(e a -6i + e b -ej) 

Yx ^ Y jx 4>% 

x^hixf 

X x <- V^Y, 

<r xy <- a xy + <j> ax (/)* y (X+X x -X y ) 
end for 
end for 



Fast MP2: 
for each x and y do 

<- G%Giy 
A*><— fi^T* 
end for 

for each a do 
for each ;' do 



"xy "xy 

end for 
end for 



4V0 
for each e do 



G XZ X Z y 

-j e °+n e niY 



e* L ai ± xy 



Yxy ^ Vxz^zy 

Xy <— Vy z J zz 

for each x and y do 

X- (Yy X -Xy)Gjy 



end for 
end for 



„a0 
Axy 



Fast BRPA: 

A^(new) <- 

for each e do 

4/ 



X xy <- fyi4&(61d)(l-^)/(a; e -w / + ^) 



^XZ^Zy 



x, v ^A^(oid)y z _ y 



^-xy 

A^(new)^A^(new) + y«B 

^xy «~ ^fi^ 

A{ y (new) A{ y (new) + Vt e Y xy j (w e - co f + S e f Q e ) 



Y xy <- O e V e ,B; 



A^,(new) <- A^(new)+X xz 7 zy 



end for 



for each e do 

4/ 



«- n f A J xy (l-5 e f )/(uj e -u f + S e f ) 



Y X y ^T- X xz V Z y 
X X y A xz V z y 



B xy <-B% 



ixy 



XxAfy 



-Y A e0 

J xz A zy 



B f xy 4- B f xy - fl e Y xy ( 1 - <f«)/(w e - W/ + $«) 



fixv^^xv 

end for 



■X xz Y zy 



for each z and e do 
for each x do 

Xy <— V^XyGxy 

^y ^ ^ywX\v 

r^xy *r- JS. X y ~T~Ay I y 

end for 

for each y do 

X x i X x G xy 
Y x 4 VxwXw 
Y x ^V zy (Y x + J$) 

end for 
end for 



'xy 



Ax 2 Cr z 



aO 



xzA z> . 
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for each a do 

X xy <s— h xz Gj y 
for each / do 
j' A j— r'- y 
end for 
end for 

K^O 
for each e do 

x xy ^Alv zv n e 

for each / do 

K X y 4- A^J, + Si^X^y 

end for 



X' 



/ 



Xy 

Yxy Y X y — X x 

for each a do 

Xjjj, <— K X y + n^JQy 

end for 
end for 

for each z do 
for each x do 



4 GxwGwzVwv 
-fb v bj 



x:'- av.Vv 

for each y do 



end for 
end for 
end for 



Kly+nijX'i 



for each z do 
for each y do 

x: h - x:'x: j 

for each x do 

end for 
end for 
end for 

II. H : DISSOCIATION 



For the model in Eq. ( 22 1, we choose reference orbitals 



0: 



i j a b 
10 10 
10-1 
10-10 
10 1 



(1) 



where {i,j} label the two occupied orbitals and {a,b} label 
the two virtual orbitals. The nonzero matrix elements are 



h\ = hj = ^-t, h a a = h b b = ^ + tl (2) 
C = ^ = -V>i(f/-V), 



where the last case applies to all Vpq matrix elements. 
The total energy and mean-field BCCD equations are 

E=E +2([i-t) + V+±(U-V)(l-T), (3) 
uj = 2t + U-(U-V)(\-T), 



where T = and u> is the orbital energy difference between 
occupied and virtual orbitals. These are also the complete HF 
equations if T = 0. In BCCD theory, T is the solution to 



= 2(uj-U)T-UU-V)(1-4T + 3T 2 ). 



The RPA+SOSEX approximation of this condition is 



(4) 



Q = 2loT-\(U-V)(\-AT + AT 2 ), (5) 



which reduces the 4-by-4 matrix RPA Riccati equation on 7y 
to a direct condition on T by using matrix symmetries. Further 
solution of these equations is straightforward algebra. 
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